Random matrix model for QCD3 staggered fermions 
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We show that the lowest part of the eigenvalue density of the staggered fermion operator in 
lattice QCD3 at small lattice coupling constant /? has exactly the same shape as in QCD4. This 
observation is quite surprising, since universal properties of the QCD3 Dirac operator are expected 
to be described by a non-chiral matrix model. We show that this effect is related to the specific 
nature of the staggered fermion discretization and that the eigenvalue density evolves towards the 
non-chiral random matrix prediction when /3 is increased and the continuum limit is approached. 
We propose a two-matrix model with one free parameter which interpolates between the two limits 
and very well mimics the pattern of evolution with /3 of the eigenvalue density of the staggered 
■ fermion operator in QCD3. 
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Introduction 



One can argue, referring to the universality [1], that the low energy properties of the QCD4 Dirac operator in the 
e- regime are described by the chiral random matrix model [2]. In particular, the microscopic eigenvalue density of the 
Dirac operator, which is obtained from the eigenvalue density 

p(A) = (^5(A-A,)) (1) 



by blowing it up at A = 0, is expected to have the same universal shape as in the random matrix. More precisely, the 
^ ' microscopic eigenvalue density is defined as 

o . 

\^ . where S — 7rp(0) and TV is the number of eigenvalues' pairs of the underlying discretized (regularized) Dirac operator, 
N is related to the physical volume V ~ Na^ where a is an UV-cut-ofF. The microscopic density p*(A) can be 
, determined analytically using the random matrix model. For the trivial topological sector and for Nf — flavors it 
reads [2] 

: P4*(A) = ^(Jo^(A) + J^(A)) . (3) 

b 

^ This eigenvalue density has been compared to numerical data from quenched Monte-Carlo simulations of lattice QCD 
with staggered fermions [3] . One could indeed see a very good agreement between the random matrix prediction and 
lattice data [4] . It has later been shown that with an improved staggered action also the non trivial topological sectors 
are well described [5]. 

If one applies the same arguments to QCD3 one is led to a non-chiral random matrix model and to the following 
microscopic density [6] 

P3*(A) = - . (4) 

TT 

Also this prediction can be easily tested numerically by comparing it to lattice QCD3 data. Actually some tests have 
been done in the past but they were not conclusive [7] . In this work we repeat quenched simulations of lattice QCD3 
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FIG. 1: Comparison of the lowest part of the eigenvalue distribution of the staggered fermion operator obtained in Monte-Carlo 

simulations of quenched QCD3 with the SU{S) Wilson action for /? = 6 on a 14^ lattice to the chiral random matrix prediction 
(blue dotted line) and the corresponding non-chiral random matrix prediction (green solid line). The data fits perfectly to the 
chiral random matrix density. 

to determine the eigenvalue density of the staggered fermion operator. The spectrum of the lowest positive eigenvalues 
of the staggered fermion operator for quenched QCD3 with (3 = 6.0 is shown in figure 1. As one can see the result is 
quite surprising. The QCD3 data matches the chiral prediction (3) and not the non-chiral one (4) anticipated in this 
case. This finding raises the following questions for discussion. Why a three dimensional theory which has no chiral 
symmetry shows a pattern characteristic for a chiral theory? Is the effect related to the regularization scheme? If yes, 
will it disappear when the continuum limit is approached? In QCD4 where the chiral symmetry plays an important 
role, the pattern will stay intact when /3 is increased and the continuum limit is approached while in QCD3 the shape 
of the lowest part of the spectrum should gradually evolve from (3) to (4). These questions will be addressed in the 
remaining part of the paper. The paper is organized as follows. In the next section we recall the matrix models used to 
describe the universal properties of the lowest part of the Dirac operator spectrum in four and three dimensions. We 
emphasize the difference between them. Then we discuss staggered fermions and introduce a new matrix model with 
one free parameter to imitate the evolution with f3 of the staggered fermion operator spectra in QCD3. It describes 
a gradual disappearance of a hard edge at the origin of the spectrum when P goes to infinity. We compare the low 
energy part of the eigenvalue distribution of quenched QCD3 and of the random matrix model. We conclude the 
paper with a short summary. 

Random matrix model 

The relation of the low energy spectrum of the Dirac operator in four dimensions to the chiral random matrix model 
has been extensively discussed in the literature so here we will restrict ourselves to a minimal presentation which is 
concentrated on the difference between the chiral and the non-chiral case. The interested reader is referred to reviews 
[8, 9] and the references therein. In four dimensions the Euclidean Dirac operator assumes in the chiral representation 
the following form 



(5) 
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where 

3 

w = mi + z^WfeCTfe , (6) 
fe=i 

1 is a 2 X 2 unit matrix, cjfc, k = 1,2,3 are Pauli matrices. The operators = —18^ + n = 1,2,3,4 are 

hermitian. In addition to the space-time indices they contain also SU{3) group indices which are suppressed in our 
notation. 

The corresponding random matrix model is obtained by replacing W hy an N x N complex random matrix D that 
inherits the block structure and basic symmetries of the operator V 

The simplest candidate for the matrix W is a Gaussian random matrix generated with the probability measure 

dfi^iW) = dWe-^^' ^^'^ , (8) 

where dW = dReWjj dlmWjj . This choice may look arbitrary at first sight but luckily it is known that the 
microscopic properties of the eigenvalue density of random matrices exhibit a large degree of universality which means 
that the microscopic density does not change for a large class of probability measures as long as they do not change 
the symmetry or the nature of the noise of the random matrix or introduce strong long range correlations. For 
example one can rigorously prove that the microscopic eigenvalue density of any random matrix generated with a 
probability measure dW exp—NTrV{W'^W) for a polynomial potential V is exactly the same as for the Gaussian 
measure (8). The microscopic properties can be however modified by adding non-analytic terms to the potential like 
for instance logarithmic ones. Actually such logarithmic terms come naturally into play in QCD if one integrates out 
fermionic degrees of freedom. Here we will however restrict ourselves only to the quenched approximation where such 
terms arc absent. In this case the Gaussian measure is the simplest and the best candidate. The behavior of the 
lowest part of the spectrum can be also modified when one replaces a square matrix by a rectangular N-^. x N- 
one, where + ~ 2N. Such a modification introduces ~ iV+ — 7V_ right-handed zero modes to the matrix 
D, if > N_, or v = N_ — N_^_ left-handed ones, if N_ > N_^_, imitating different topological sectors. In the 
trivial topological sector, which we consider here, W can be viewed as a square Gaussian complex matrix (8). The 
microscopic eigenvalue distribution of the matrix D (7) can be calculated analytically. The calculations [2] give the 
expression mentioned above in eq. (3). 

Consider now three dimensions. In this case there is no chiral symmetry as we will see below but one can define 
a two fcrmion family model which imitates this symmetry. The spinors have two components in the fundamental 
representation. There are two independent nonequivalent sets of gamma matrices jk = o'k and = — cxfe, where as 
before crj, are Pauli matrices. One can associate a fermion family with each representation and consider a common 
theory for the two families. A Lagrangian for such a theory can be concisely written in terms of four component 
spinors [10]. It takes a very similar form to the four dimensional Lagrangian. In particular it has a U(2) symmetry, 
which is broken by a mass term to a U{1) (8) U{1) symmetry. We thus have a situation similar to the chiral symmetry 
in four dimensions. The main difference is that the term with 74 is absent in the Lagrangian. The fermionic operator 
takes the form 



W = (9) 



fc=i 



So in three dimension the operator W is anti-hermitian W = — W^. Thus if one constructs a matrix model one should 
substitute a complex matrix W by an anti-hermitian one. Let us write W = iB 



^3=( ,H^t =U ) (10) 



where B = is from a GUE ensemble with the standard measure 

df,s{B) = dBe-^^^' , (11) 

and where dB = Ylu dBa Y[i<j dKcBijdlmBij is a flat measure in the set of hermitian matrices. As for D4 (7) the 
eigenvalues of the matrix D3 are purely imaginary. They also come in pairs ±iX. The microscopic eigenvalue density 
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is however different. It is given by the expression (4), so it has neither a dip nor a wavy structure. This shape just 
resuhs from blowing up the central region of the Wigner semicircle eigenvalue distribution of the matrix B which is 
flat and has no hard edge at A = 0. One expects that the fermionic operator of QCD3 should reproduce this structure 
in the continuum limit. So we have to understand the eigenvalue density observed for staggered fermions for quenched 
QCD3 shown in figure 1. 

Staggered fermions 

Staggered fermions discretization and its relation to continuum physics in even dimensions was discussed in [11] and 
in odd dimensions in [12]. Let us shortly discuss how to calculate of the spectrum of the staggered fermion operator 
in quenched simulations of lattice QCD. In our simulations we used the standard Wilson action to generate gauge 
fields 




where U^{x) is a SU{3) matrix - a dynamical variable associated with a link going between two neighboring sites x 
and X + II. The first sum runs over all lattice sites x, the second one over all pairs of directions I < < v < D, so 
Sw collects contributions from all elementary plaquettes. In the quenched approximation one generates gauge field 
configurations with the probability measure JJ dUe~^'^ which is independent of fermionic degrees of freedom. One 
uses these configurations to compute quantum averages. In other words the influence of fermions on gauge fields is 
neglected in this approximation. 

The staggered fermion operator is defined as [3] 

d 

D% = E^'^(^) K'(^x-+>^ - (u^y^i^ - m)<5,,.+m) (13) 

IJ.=0 

where ij are indices of the SU{3) matrix and r]fj,{x) = - for = 1, . . . , D. One can also define a parity 

operator e{x) = ■+^d ^ jf the lattice has even number of sites in each direction, that divides sites into odd and 

even ones and introduces a chessboard structure. Each odd site has even neighbors and vice versa, so the staggered 
fermion matrix can be decomposed into blocks 

which are mutually anti-conjugated D^o = —D\^ by construction (13). Writing Doe = iW, one obtains a matrix of 
the form (10). Figure 1 tells us that the Dgo = iW behaves for (3 = 6.0 as if W belonged to the universality class 
of Gaussian complex matrices. Since the even-odd split Df-o = is built into the discretization scheme one can 

expect that the related microscopic universality is robust unless something dramatic happens that changes the basic 
symmetry or the nature of the randomness of the matrix. An example of such a mechanism may be a freeze-out of 
some degrees of freedom of the matrix. 

In fact, it has been shown that in the continuum limit the staggered fermions in three dimensions corresponds to 
four two component spinors with the action 

SF = ^[u{ak®I)D^u-d{(Jh'»I)D^,d\ (15) 

where / is a two by two unit matrix. This can be equivalently written as two copies of the four component fermions 
defined above. Thus we expect that in the continuum limit, that is for /3 — > 00, the matrix W becomes anti-hermitian, 
or equivalently that the fluctuations of hermitian degrees of freedom are suppressed in this limit. In the next section 
we propose a random matrix model describing such a gradual suppression. This model has one parameter which 
interpolates between the regime where is a generic random complex matrix (with both the hermitian and anti- 
hermitian sectors) and the regime where it is an anti-hermitian one W = —W^. As it is shown below the model 
captures the pattern of evolution with /3 of the microscopic density observed in the numerical QCD3 simulations. 
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Freeze-out random matrix model 



The idea is to consider a random complex matrix W which is a linear combination W = xA+iyB of two independent 
identically distributed Gaussian hermitian matrices A = and B = with real coefficient x, y. The matrices are 
generated with the probability measure 

dti{A, B) = d^dSe-^^^'e-"^^^' , (16) 

where dA and dB are flat measures for hermitian matrices. To avoid redundancy with the width parameter E it 
is convenient to restrict x,y to x"^ + = 1 or equivalently to parametrize x = cos(q!), y = sin(a). We obtain a 
one-parameter family of matrices 

„ _ / iW'<\_f icos{a)A-sm{a)B\ , . 

"~\iW J ~ \ icos{a)A + sm{a)B J ' '' '' 

with a mixing parameter a which interpolates between (7) and (10). For a = 7r/4 this matrix is equivalent to (7) 
where W = {A + iB)/^/2 (8). The integration measure for W can be derived by changing variables in (16) which 
gives: 

dii{W) ~ dWe-^^'^^^' (18) 

This is a standard Girko-Ginibre ensemble [13, 14]. Eigenvalues of W are uniformly distributed in a disk of radius 
l/S in complex plane centered around the origin. For arbitrary a it becomes an elliptic ensemble with a measure 

d^,{W) ^ We-(^^(^^'-*(^^+^'^')) (19) 

where r = cos(2a) [15]. Eigenvalues are now uniformly distributed in an elliptic disc with semi-axes of relative length 
(1 — t)/(1 + r). When a approaches 7r/2, hermitian degrees of freedom are gradually suppressed. For a = 7r/2, W 
becomes purely anti-hermitian (10) and the ellipse reduces to a cut on the real axis. The evolution of the eigenvalue 
spectrum of the matrix Da for a from 7r/4 to 7r/2 (or equivalently r from to —1) smoothly interpolates between the 
limiting eigenvalue densities /34*(A) and pz*{^) when a. 

It is instructive to begin the analysis of the spectrum of the matrix Da by considering the case A'' = 1. Although it 
is a very simplified version of the full model it allows one to understand the nature of the hard edge of the spectrum 
and the mechanism of its disappearance for a — > 7r/2. For = 1 the matrices A, B reduce to real numbers a, b which 
are independent, identically distributed normal random variables with the variance 1/S^. The matrix (17) has only 
one pair of eigenvalues 

A± = ±ir = ±i^a? cos'^{a)+h^ sm^{a) . (20) 

We are now interested in finding the probability distribution of the random variable r which is constructed from a 
and b. This quantity has a geometrical interpretation. The vector r = (cos(a)a, sin(a)6) can be viewed as a random 
vector on a plane whose components are independent Gaussian variables with the standard deviations cos(a) 
and yj~^ sin a. The length of this vectors is equal to r (20). The probability distribution for this vector is an elliptic 
Gaussian (lines of equal probability form ellipses) so it is easy to derive the corresponding distribution of length r of 
the vector f. For a = 7r/2 the distribution reduces to a one dimensional distribution concentrated on the 6-axis. In 
this case r = \b\ and thus the probability distribution of r is a half-Guassian distribution p(r) = 2S/\/27re^^ For 
a = Tr/4 the distribution is a spherically symmetric Gaussian one and the distribution of r can be easily calculated 
by changing to the polar coordinates and integrating out the angle. It yields p{r) = 2Y?re~^ . For a in-between 
the distribution can be expressed by the following integral 

^ 2tt Jq cos^(a) cos^(0) -|- sin^(a) sin^(^) ^\ 2(cos^(a) cos^(^) -|- sin^(a) sin^(^)) / ^ ^ 

which interpolates between the two limits. The evolution of the probability distribution with a is shown in figure 2. 
For small r the integral behaves as Pa{r) ~ 2E^r/ sin(2a), so it goes linearly to zero with a slope 2S^/ sin(2a). The 
slope becomes infinite for a — )■ 7r/2. In this limit also the position of the maximum goes to zero and the dip gets 
narrower. It disappears completely only for a = 7r/2. To summarize, we see that as long as a 7^ tt/2 there is a 
repulsion from the hard edge at the origin. It can be attributed to the two-dimensional nature of the underlying 
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FIG. 2: The distribution pi (A) of the lowest eigenvalue for the systems considered in this article. (Left) The distribution (21) 
Pa{r) for cot(a) = 1.0, 0.4, 0.2, 0.1, 0.5, 0.0. We set S = 1. The maximum moves from right to left when a is decreased. The 
function Pa{r) goes linearly to zero Pa(r) ~ 2r/sin(2a) for r ^ with a coefficient 2/sin(2a). The dip disappears only for 
sin(2a) — 0. (Middle) The distribution of the lowest eigenvalue obtained numerically for the matrix Da in the matrix model 
for A'^ = 20 and for cot(a) = 0.0, 0.02, 0.05, 0.10, 1.0. As before the maximum moves from right to left when a is increased from 
7r/4 to 7r/2. The linear part of the slope at zero has a coefficient which increases when a approaches -k/2. (Right) QCD3 for 
/? = 12 and 30 on 24^ lattices. The eigenvalues are rescaled according to (2). 

distribution of the vector {cos{a)a, sin(a')6) in the plane. The hard edge and the repulsion completely disappears only 
for sin(2a) = where the distribution becomes one-dimensional. We expect the same mechanism of repulsion for the 
smallest eigenvalue of the matrix (17) also for A'' > 1. As long as both hermitian and anti-hermitian degrees of 
freedom are active the lowest eigenvalues will be repelled from zero. Indeed we sec in figure 2 in the middle that the 
evolution of the probability density for the lowest eigenvalue of the matrix Da for N = 20 very much resembles that 
shown in the left panel. For large a the distribution has a shape of the type re~'^ while for small ones it is more 
like a semi-Gaussian shape with an additional narrow peak on top of it. This additional peak gets narrower when 
a — ^ 7r/2 and it moves towards zero. The same pattern is observed also for larger TV. Finally, in the right chart we 
show corresponding plots in QCD3 data. 

The repulsion has also an influence on the second smallest eigenvalue, the third smallest one etc. For a — > 7r/2 the 
repulsion range becomes smaller and correspondingly tlic influence on higher eigenvalues becomes smaller too. The 
resulting microscopic eigenvalue distribution will be discussed in the next section where it will be compared to the 
corresponding histograms for quenched QCD3 for different values of /3. 

Comparison of QCD3 to the random matrix model 

In figure 3 we show histograms of eigenvalue distribution for QCD3 for different values oi /S — 6.0; 12.0; 18.0 for 24^ 
lattice (Ihs) and for the matrix model (5) for different values of cot(a) = 1.0; 0.1; 0.05 and N = 100 (rhs). When /3 and 
a decrease the wavy structure gradually disappears in a very similar way in both cases as one can see by comparing 
the figures. One should notice an anomaly in the plot in the bottom-left panel for /3 = 18.0 where the main peak 
seems to be slightly wider than the corresponding peaks in other panels. The broadening is related to an effect of 
eigenvalue pairing which appears when /3 is increased. The pairing means that the distance between the first and 
the second eigenvalue, the third and the fourth one, etc decreases when /? increases. For /? = 18.0 the distribution of 
the first one and the second one are so close to each other, as shown in figure 4, that they contribute to one broader 
peak in the eigenvalue density function seen in the bottom-left plot in figure 3. For /3 ^ 00 the distance between 
neighboring eigenvalues tends to zero giving degenerate pairs of eigenvalues in the limit as one can see in the right 
chart in figure 4. This is in agreement with what one expects from the continuum theory (15) for which the Dirac 
operator indeed has doubly degenerate eigenvalues. 

Summary and discussion 

We showed that a simple two-matrix model reproduces basic characteristic features observed in the eigenvalue 

spectra of the Dirac operator in staggered QCD3 data. The mixing parameter r = cos 2a interpolates between the 
regime where block matrices W (19) in the Dirac matrix (17) are complex to the regime where they are anti-hermitian. 




FIG. 3: Left: eigenvalue distribution for QCD3 for /3 = 6.0; 12.0; 18.0 for lattices 14^ , 24^, 24^ respectively. Right: eigenvalue 
distribution for Da for the matrix model for cot(a) = 1.0; 0.1 and 0.05 for 50 x 50 matrices. 



In the former case the microscopic density is given by (3) while in the latter one by (4). We observe a similar evolution 
in the QCD3 data. The shape of the probability distribution of the lowest eigenvalue evolves in both the cases in 
a similar way too. The nature of finite size effects manifesting as the appearance of a small narrow peak on top of 
the distribution as well as the gradual disappearance of the hard core close to the origin is similar in both the cases 
too. What remains to show is how the microscopic density precisely interpolates between the two limiting cases (3) 
for Q: = 7r/4 (t — 0) and (4) for a = n/2 (r = —1). We expect that a continuous evolution of the microscopic 
density can observed in a parameter a = ^/N{tt/2 — q) in a double scaling limit, that is for A'' — >■ 00, a ^ 7r/2 and 
for finite a. This limit corresponds to a very asymmetric mixing of hermitian and anti-hermitian degrees of freedom 
W ss {a/^fN)A + iB. The width of fluctuations of hermitian degrees of freedom decreases with the system size as 
1/VN and the aspect ratio r of the elliptic ensemble (19) approaches minus one as t = — 1 + Aa^/N. The evolution 
of the shape with a can be found using for instance the supersymmetric method [16]. Having done that one can try 
to relate the parameter S to ^ in QCD3. We leave this issue for future investigations. 
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FIG. 4: (Left) Histograms of individual eigenvalues for QCD3 with /3 = 18.0 on 24^ lattice. (Right) Averages of the first six 
eigenvalues in QCD3 for N = 24^ as a function of /3. All eigenvalues are normalized to the smallest eigenvalue for given /3. 
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